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We describe a fluctuating volume-current formulation of electromagnetic fluctuations that extends our recent 
work on heat exchange and Casimir interactions between arbitrarily shaped homogeneous bodies [Phys. Rev. 

B. 88, 054305] to situations involving incandescence and luminescence problems, including thermal radiation, 
heat transfer, Casimir forces, spontaneous emission, fluorescence, and Raman scattering, in inhomogeneous 
media. Unlike previous scattering formulations based on field and/or surface unknowns, our work exploits 
powerful techniques from the volume-integral equation (VIE) method, in which electromagnetic scattering is 
described in terms of volumetric, current unknowns throughout the bodies. The resulting trace formulas (boxed 
equations) involve products of well-studied VIE matrices and describe power and momentum transfer between 
objects with spatially varying material properties and fluctuation characteristics. We demonstrate that thanks 
to the low-rank properties of the associated matrices, these formulas are susceptible to fast-trace computations 
based on iterative methods, making practical calculations tractable. We apply our techniques to study thermal 
radiation, heat transfer, and fluorescence in complicated geometries, checking our method against established 
techniques best suited for homogeneous bodies as well as applying it to obtain predictions of radiation from 
complex bodies with spatially varying permittivities and/or temperature profiles. 


I. INTRODUCTION 

Quantum and thermal fluctuations of charges give rise to a 
wide range of electromagnetic phenomena; these include lu¬ 
minescence from active media, e.g. fluorescence and spon¬ 
taneous emission, the finite linewidth of lasers near thresh¬ 
old,'^’^ thermal radiation and heat transfer from hot objects,®^''^ 
and dispersive interactions (Casimir forces) between nearby 
surfaces.Fluctuation-driven effects are not only responsi¬ 
ble for many naturally occurring processes but are also poised 
to take an increasingly active role in emerging nanotech¬ 
nologies,'^’*^ spurring interest in the study and engineering 
of complex shapes that could dramatically alter their behav¬ 
ior.*'*’^' Although rooted in similar principles, the physical 
mechanisms behind each of these processes vary considerably, 
leading to theoretical descriptions that differ both in their for¬ 
mulation and implementation. Ultimately, however, all such 
calculations reduce to a series of classical scattering prob- 
lems^^’^-* that until recently remained largely specialized to 
situations involving simple, high-symmetry geometries, e.g. 
planar and spherical objects. 

In this manuscript, we present a framework for the general- 
purpose calculation of many different incandescence and lu¬ 
minescence processes, including fluorescence, spontaneous 
emission, thermal radiation, heat transfer, and Casimir forces 
in arbitrary geometries. In particular, we derive a fluctuating 
volume-current (FVC) formulation of electromagnetic fluctu¬ 
ations that exploits techniques from the volume-integral equa¬ 
tion (VIE) formulation of electromagnetic scattering^'*’^^ and 
which expands the range and validity of current methods to 
situations involving inhomogeneous media. Although FVC 
is similar in spirit to our previous fluctuating surface-current 
(FSC) methods,unlike FSC our new approach is not lim¬ 
ited to piecewise-homogeneous objects. Here, the unknowns 


are volume currents within objects rather than surface cur¬ 
rents as in FSC, and can therefore easily handle more com¬ 
plex structures, including inhomogeneous bodies with tem¬ 
perature gradients or spatially varying permittivities. In con¬ 
trast to recently developed scattering-matrix methods,the 
FVC and FSC methods do not require a separate basis of in¬ 
coming/outgoing wave solutions to be selected (a potentially 
difficult task in geometries involving interleaved objects or 
complex structures favoring nonuniform spatial resolution), 
although VIE can be used to compute the scattering matrix 
if desired. We show that regardless of which quantity is com¬ 
puted, the final expressions for power and momentum trans¬ 
fer are based on simple trace formulas involving well-studied 
VIE and current-current correlation matrices that encode the 
spectral properties of fluctuating sources. We And that while 
the number of VIE unknowns is large compared to scattering 
or ESC formulations, the associated VIE matrices admit low- 
rank approximations that turn out to significantly reduce the 
complexity of trace evaluations, making practical calculations 
tractable. We validate the EVC method by checking its predic¬ 
tions against known solutions for homogeneous objects and 
then apply it to calculate thermal radiation, heat transfer, and 
fluorescence from compact objects (spheres, ellipsoids, and 
cubes) with spatially varying permittivities and temperature 
gradients. The same trace formulas can be readily adapted to 
obtain the angular distribution of far-fleld radiation, which we 
illustrate by providing new predictions of directional emission 
from inhomogeneous objects, the As explained below, while 
VIE methods can be applied to arbitrary geometries, they are 
particularly advantageous in situations where object sizes are 
on the order of (or smaller) than the relevant wavelengths, 
providing a useful complement to well-established techniques 
better suited for the study of arbitrary geometries with length- 
scales that are large or small compared to the relevant electro- 
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magnetic wavelengths, e.g. proximity approximations.'®"^^ 

Electromagnetic fluctuation phenomena can be roughly di¬ 
vided into two categories: incandescence and luminescence 
problems. Incandescence refers to electromagnetic radiation 
from objects generated by the quantum and thermal motion 
of charged particles in matter, whereas luminescence refers to 
incoherent emission of light from non-thermal sources. The 
oldest and most well-studied manifestation of incandescence 
is the familiar glow of objects—thermal radiation—that oc¬ 
curs when an object is heated above the temperature of its 
surrounding environment."'^’"'^ Although Planck’s law was not 
more than a century ago at the center of vigorous contro¬ 
versy which helped establish the foundations of quantum me¬ 
chanics,"'"' much of our recent interest in this phenomenon 
spawns from its profound impact on energy and related nan¬ 
otechnologies. Interest in complex designs is also fueled by 
our increasing ability to engineer selective and even dynami¬ 
cally tunable emitters and detectors at wavelengths for which 
there is currently a lack of coherent sources,in addi¬ 
tion to solar-energy harvesting applications.In addition 
to radiation, fluctuations can also mediate heat exchange®’*’^^ 
and interactions®’'^’^'’^®’^® (known as Casimir forces) between 
objects—unlike heat exchange, Casimir interactions persist 
even at equilibrium and are known to arise primarily due to 
contributions of quantum rather than finite-temperature fluctu¬ 
ations. One fundamental distinction between “near-held” ef¬ 
fects (between objects at wavelength-scale separations or less) 
and the more familiar “far-held” phenomena (separations ^ 
wavelength) is that the former can be signihcantly enhanced 
by the contributions of evanescent waves,growing in 
a power-law fashion with decreasing object separations. As a 
result, the heat transfer between real materials can exceed the 
predictions of the Planck blackbody law by orders of magni¬ 
tude'^ and quantum forces can even reach atmospheric pres¬ 
sures at nanometric lengthscales,^' motivating interest in com¬ 
plex designs that can be tailored for various applications, in¬ 
cluding thermophotovoltaic energy conversion,®'’^®-' nanoscale 
cooling,®"'’®® and MEMS design.®'^®* 

Until very recently, however, calculations and experiments 
remained focused on planar structures and simple approx¬ 
imations thereof.®“'"'’®® Since all such thermal effects arise 
due to the presence of huctuating current sources, from the 
perspective of calculations their descriptions reduce to a se¬ 
ries of classical scattering calculations involving helds due 
to currents,'"'’^® the spectral characteristics of which are re¬ 
lated to the underlying physical means of excitations. In the 
case of incandescence, they are determined by the thermal 
and dissipative properties of materials via the well-known 
fluctuation-dissipation theorem (EDT).®"’®' Naively, this in¬ 
volves repeated calculations of electromagnetic Green’s func¬ 
tions throughout the bodies, which can prove prohibitive for 
complex objects where the latter must be computed numeri¬ 
cally, especially due to the broad bandwidth associated with 
thermal fluctuations, but it turns out that more sophisticated 
formulations exist.These include time- and frequency- 
domain methods where the power transfer or force on an ob¬ 
ject is obtained via integrals of the flux or Maxwell stress ten¬ 
sor, or equivalently electromagnetic Green’s functions, along 


some arbitrary surface enclosing the body.®®’®^^®* Recent tech¬ 
niques forgo surface integrations altogether in favor of unfa¬ 
miliar but more efficient expressions involving traces of ei¬ 
ther scattering®'’®®’®"''®®^®®'®®’*'’ or boundary-element®®’®®’*' ma¬ 
trices. Regardless of the choice of unknowns, in practical im¬ 
plementations the latter are expanded in terms of either delo¬ 
calized spectral bases (e.g. Eourier or Mie series) best suited 
for high-symmetry geometries, or geometry-agnostic local¬ 
ized bases (piecewise polynomial “element” functions) de¬ 
fined on meshes or grids and applicable to arbitrary objects.®® 
While there has been much progress so far, these methods 
have yet to be generalized to handle structures with temper¬ 
ature gradients or varying permittivities. 

Temperature gradients can arise for instance due to the in¬ 
terplay of phonon and photon transport,*® *® such as in hetero¬ 
geneous structures with disparate thermal conductivities, in¬ 
cluding chalcogenide/metal interfaces*"' *® or quartz-platium- 
polymer structures,*® or in graphene-based devices.*® Tem¬ 
perature gradients have also been observed in atomic force 
microscopes**’*® and nanowires,®'’ as well as in situations in¬ 
volving irradiated particles immersed in fluids,®’“"’'’ magnetic 
nanocontacts,"” or microcavities subject to strong photother- 
mal effects."’® Material inhomogeneities also arise in mi¬ 
crocavity lasers stemming from nonlinear effects."’® Surpris¬ 
ingly, there are only a handful of calculations involving non- 
isothermal particles, including calculation of radiation from 
atomic gases in shock-layer structures with linear tempera¬ 
ture gradients"’"' or calculations of large-radii spheres based 
on Mie series or related semi-analytical expansions."’®'"’® As 
we show in a separate publication, temperature gradients in 
inhomogeneous bodies can lead to a number of interesting ef¬ 
fects, including highly directional thermal emission."’® 

Luminescence, like incandescence, involves incoherent 
emission of light due to quantum and thermal fluctuations 
of charges, but differs in that excitations are driven by co¬ 
herent rather than thermal sources. Examples include spon¬ 
taneous emission, Raman scattering, and fluorescence from 
active media externally pumped by coherent light. ®’"’*’"’® Al¬ 
though the spectral properties of fluctuating currents depend 
on complicated and often nonlinear light-matter interactions, 
the resulting radiation is incoherent and can be modeled by 
exploiting scattering techniques similar to those employed in 
incandescence problems."’® There are however many impor¬ 
tant differences between these two classes of problems. Eor 
instance, the luminescence spectrum of many emitters is rel¬ 
atively narrow (involving wavelengths close to material reso¬ 
nances) and this has implications for calculations which favor 
frequency as opposed to time-domain techniques (the latter 
being better suited for broad-bandwidth processes). Eurther- 
more, while many thermal radiation problems involve objects 
with uniform temperature distributions, the properties of cur¬ 
rent fluctuations excited by external pumps depend sensitively 
on the inputs and can change dramatically and continuously 
throughout the bodies, which is problematic for SIE/ESC for¬ 
mulations based on piecewise homogeneity. Such a situation 
arises for instance in the fluorescence from objects with fea¬ 
tures ^ incident wavelengths, where resonant absorption can 
lead to significant spatial variations in the amplitudes of the 
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fluctuating currents.^ 

Until recently, the fluorescence or Raman emission pat¬ 
tern of small particles was obtained by analytical methods 
based on Mie series or related basis expansions.*'*^’*" More 
recent techniques for studying luminescence from arbitrarily 
shaped particles instead rely on numerical techniques,**^ most 
commonly time-domain methods,**^ *** and include stud¬ 
ies of bowtie antennas,**® nanostars,*^** conical tips,*^*“*^^ 
dimers,*^"* and thin fllms.*^^ Frequency domain methods in¬ 
clude finite-element,*^®'*^® boundary-element,*®* and discrete 
dipole approximation (DDA)*®® *®® methods. These tools 
have been exploited for instance to demonstrate that both 
shape and material degrees of freedom can be used to tailor 
particle emission, making it possible to enhance fluorescence 
and Raman processes®'****’**’® as well as obtain unusual angular 
emission patterns;*®® *®® even more recently, there has been 
interest in studying effects related to active (non-Hermitian) 
systems.*®®“*®® In most cases (with a few exceptions**®), the 
total radiated power in a given direction is computed by di¬ 
rectly summing the contribution of individual emitters inside 
the objects, requiring repeated evaluation of Green’s functions 
over both volumes and surfaces. In addition, many calcula¬ 
tions rely on approximations in which the effect of the in¬ 
cident drive is either approximated or entirely neglected*"*** 
or where only the radiation from a partial set of emitters in¬ 
side the objects is obtained.*"** Our FVC-VIE approach not 
only removes limitations associated with such approxima¬ 
tions by fully accounting for both the emission and excitation- 
dependent properties of all fluctuating sources, but introduces 
new trace-formulas that offer compactness, simplicity and a 
unified framework for computing a wide range of fluctuation 
phenomena, allowing techniques and ideas from one area to 
be more easily applied to another. 

A technique that in principle shares many similarities with 
the VIE method is the so-called discrete-dipole approximation 
(DDA),*"*® which models objects as finite arrays of polarizable 
dipoles whose response and interactions due to incident elec¬ 
tromagnetic fields can be obtained via the solution of a cor¬ 
responding integral equation. *"*® DDA has been recently em¬ 
ployed and suggested as an efficient approach for computing 
radiative heat transfer*"*"* as well as fluorescence®’*"*® from ar¬ 
bitrary geometries, but unfortunately suffers from a number of 
important limitations. Technically, DDA belongs to the gen¬ 
eral class of volume integral equations traditionally solved nu¬ 
merically via the method of weighted residuals*"*® (or method 
of moments as it is conventionally known when applied to 
computational electromagnetics*"*®), by which integral equa¬ 
tions are converted into a solvable and finite set of linear sys¬ 
tems of equations. Specifically, system unknowns (fields or 
equivalent currents) are approximated by expanding them in 
a finite set of basis functions, often determined by discretiza¬ 
tions of objects into meshes or grids, and then forcing the re¬ 
sulting semi-discrete equations to be equal in a weak sense, 
i.e. by integrating them against a set of testing functions.*"*® 
The actual choice and combination of basis and testing func¬ 
tions gives rise to a plethora of practical variants.*"*® 

DDA can be considered to be a particular implementation 
of the VIE method known as a collocation method,*"** involv¬ 


ing constant or dipole basis functions and Dirac-delta distri¬ 
butions for testing, with solutions forced to be accurate only 
at a finite set of points (known as point matching).*"** How¬ 
ever, it is now known that methods of weighted residuals 
are only guaranteed to converge in norm under special cir¬ 
cumstances, the lack of which can lead to numerous conver¬ 
gence and efficiency issues.*"*® Specifically, basis functions 
must span the function space of the unknowns and testing 
functions must span the dual space of the range of the corre¬ 
sponding VIE operator.*®**’*®* DDA respects neither of these, 
and as a consequence its applicability is largely limited to sit¬ 
uations involving light scattering in structures with small in¬ 
dex contrasts and weakly polarizable media,*"*® beyond which 
it can lead to a number of severe convergence and accu¬ 
racy problems.*®® (Note that DDA also makes a number of 
other approximations that break down in geometries involv¬ 
ing wavelength-scale objects, cf. Eq. 14 in Ref. 143.) In 
contrast, our EVC formulation is based on a recently devel¬ 
oped VIE framework (dubbed JM-VIE) that is numerically 
solved by means of a Galerkin method of moments.®® JM- 
VIE exploits basis and testing functions spanning the func¬ 
tion space of internal volume currents,®® the stability and su¬ 
perior convergence of which have been demonstrated in ge¬ 
ometries involving highly inhomogeneous objects and large 
dielectric contrasts.®® While the associated JM-VIE matrix el¬ 
ements involve complicated, expensive, and highly singular 
volume-volume integrals of homogeneous Green’s functions 
integrated against pairs of basis functions, these were recently 
shown to reduce to surface-surface integrals over smoother 
kernels that can be readily handled using specialized integra¬ 
tion techniques originally developed for SIE methods.*®®’*®® 


In the following sections, we derive our EVC formulation 
of fluctuating currents and demonstrate that it can be em¬ 
ployed to study a wide class of electromagnetic fluctuation 
effects in general geometries, with no uncontrolled approx¬ 
imations except for the finite discretization (basis). We be¬ 
gin in Sec. II with a brief review of the VIE formulation of 
electromagnetic scattering, followed by derivations of formu¬ 
las involving power and momentum transfer, as well as far- 
field radiation patterns from radiating objects. The final boxed 
expressions are described via traces of products of VIE and 
current-current correlation matrices which encode the spa¬ 
tial and spectral characteristics of the fluctuating sources. In 
Sec. Ill, we show that important algebraic properties of the as¬ 
sociated VIE and correlation matrices allow efficient evalua¬ 
tion of the trace expressions; specifically, a number of the VIE 
matrices admit low-rank approximations, enabling us to ex¬ 
ploit sophisticated and fast iterative techniques for their eval¬ 
uation. Einally, in Sec. IV the EVC framework is validated 
against known results and also applied to obtain predictions in 
new geometries that currently lie outside the scope of state-of- 
the-art techniques, such as objects subject to spatially varying 
temperatures and dielectric properties. 
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component electromagnetic fields and volume currents. 



and consider the scattering problem involving incident fields 
(pine due to a (in the absence of bodies) and scattered fields 
^scat due to reflections and scattering from objects and 
sources. Defining the 6-component volume currents 


Figure 1. Schematic of a many-body geometry in which fluctuating 
current sources give rise to radiation as well as flux and momentum 
transfer between the bodies. Also illustrated are the incident field 
(j>inc due to a single dipole source a within a body Vi along with the 
induced polarization-currents ^ throughout Vi and two nearby bod¬ 
ies, V 2 and V 3 , resulting in scattered fields (j)sca.t- The characteristics 
of the dipole sources cr (fluctuation statistics) and the permittivities 
of the bodies x (material properties) both vary within each object. 


II. FVC FORMULATION 

In this section, we begin by reviewing the VIE method 
of EM scattering and apply it to derive an EVC formulation 
of fluctuation-induced phenomena in inhomogeneous media. 
Our approach relies on the JM-VIE formulation and associ¬ 
ated Galerkin method of moments presented in Ref. 25, also 
briefly discussed. As noted above, a strategy based on SIE for¬ 
mulations is unavailable for modeling inhomogeneous objects 
since finding the radiation of a point source (the Green’s func¬ 
tion) in inhomogeneous media is nearly impossible with only 
surface unknowns. Matters are further complicated for fluc¬ 
tuation phenomena involving power or momentum transfer, in 
which case inhomogeneities in the properties of the fluctuat¬ 
ing sources (e.g. spatial variations throughout the bodies due 
to temperature or dielectric changes) must also be accurately 
accounted for. Starting with the recently developed power for¬ 
mulas,we derive compact trace expressions for the power 
and momentum transfer and far-held radiation pattern of com¬ 
plicated objects with inhomogeneous properties. Einally, we 
elaborate on special algebraic properties of the associated VIE 
and correlation matrices that allow fast computations of the 
matrix-trace formulas, making large and complicated calcula¬ 
tions tractable. 


A. Volume integral equations 

The derivations of VIEs often rely on the volume equiva¬ 
lence principle, which shares many similarities with—^but is 
significantly simpler and more easily derived than—the more 
well-known surface equivalence principle. Consider the 
system of arbitrarily shaped, inhomogeneous bodies described 
by the relative permittivity e and permeability p functions, 
depicted schematically in Eig. 1. Let (p and a denote 6- 


^ = (mJ = 

associated with bound polarization J{, and magnetization M{, 
currents inside the objects, described by the 6x6 susceptibil¬ 
ity tensor x (which for convenience also includes the permit¬ 
tivity and permeability of the ambient medium), it follows that 
the scattered field can be be written as a convolution of ^ with 
the homogeneous Green’s function of the ambient medium.^^ 
(Note that there is no assumption on x, which can describe 
both anisotropic and/or chiral media, changing only the form 
of the homogeneous Green’s function.^^^) In particular, the 
unknown scattered fields can be shown to be related to the 
free and bound currents, respectively, via convolutions (*) 
with the 6x6 homogeneous Green’s tensor of the ambient 
medium (typically free space) r(x, y) = r(x — y, 0), written 
explicitly in Ref. 26. This is the core idea behind the volume 
equivalence principle, which we review below. 

We begin by writing the total field <p = E * (cr -F $) via 
the volume equivalence principle^"^ in terms of the incident 
i'inc = r * cr and scattered (psca-t = T x ^ fields, or more 
explicitly: 

= y dVr(x,y) [cr(y)-FC(y)] (2) 

where it is clear that all of the scattering information (includ¬ 
ing material inhomogeneities) is “encoded” in the convolu¬ 
tion of the homogeneous Green’s function with the polariza¬ 
tion/magnetization current. Multiplying both sides of Eq. 2 
with —iuJX and using the definition of ^ in Eq. 1, one arrives 
at the following VIE for the induced currents 

C +*wx(r*0 =-ia;x(r*f7), (3) 

which can be solved to obtain ^ from the incident sources a. 
This is the so-called JM-VIE formulation of electromagnetic 
scattering in which the unknowns are induced currents rather 
than fields or field densities. Compared to other formulations 
based on field unknowns, JM-VIE exhibits superior perfor¬ 
mance in terms of accuracy and convergence, especially for 
objects with high refractive index. 

The operator equation above is customarily solved by re¬ 
ducing it to an approximate, finite-dimensional linear system. 
Let {ba} be some convenient set of N vector-valued basis 
functions. We can then approximate our unknowns ^ (and, for 
convenience below, the source currents a) in this basis: 

N N 

^(x) W ^ Xaba{x), (t(x) PS ^ Saba{x). (4) 

CX.— 1 OL—l 
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There are two main categories of basis functions that are used 
in the numerical solution of the JM-VIE above, known as 
spectral andMoMsub-domainhases. A spectral basis consists 
of non-localized Fourier-like basis functions whereas MoM 
sub-domain bases are localized functions obtained by dis¬ 
cretizing objects into meshes or grids of volumetric elements, 
e.g. cubes, tetrahedra, and hexahedra,'®' and defining func¬ 
tions by low-order polynomials with local support in one or 
a few elements. In this work, we resort to the second cate¬ 
gory and exploit piecewise constant basis functions defined in 
cubes, due to the flexibility they offer for modeling geome¬ 
tries of arbitrary shape.We note however that the proposed 
framework and the resulting matrix-trace formulas can also be 
evaluated using spectral bases as well. 

Finally, the semi-discrete equation is “tested” with another 
set of functions (called testing functions) to produce a linear 
system. In the Galerkin approach, the set of testing functions 
is the same with the one of the basis functions. The resulting 
Galerkin JM-VIE linear system reads 


are now well-established, fast algorithms to reduce the costs 
of such integral equation solvers.However, the abil¬ 
ity to exploit fast solvers in fluctuation EM problems is not 
a priori guaranteed since as we show below the final formu¬ 
las involve complicated traces of products of JM-VIE and re¬ 
lated matrices. In Sec. Ill, we describe a fast procedure for 
the computation of the proposed matrix-trace, which relies on 
a straightforward and easily implemented FFT-based fast al¬ 
gorithm presented in Ref. 25 that scales as Cl(iV log iV) for 
each matrix-vector product and requires 0{N) memory. 

Before concluding this section, we introduce some addi¬ 
tional definitions and notation. In particular, further below we 
exploit the so-called Green matrix G, defined as 

G «,/3 = (8) 

which involves interactions among basis functions mediated 
by the Green’s function. For n objects, the associated matrices 
and vectors can be conveniently written as: 


W-^x = {V-W-^)s, (5) 


where 


bp + * bp)) 

v^,p = (b^,bp) 

and a, P = I : N. Also, (,) denotes the standard inner prod¬ 
uct of functions {<j), ip) = J 4>*ip, with the * superscripts de¬ 
noting the conjugate transpose (adjoint) operation. Without 
loss of generality, we can choose the basis functions to sat¬ 
isfy an orthogonality relation, so that {be, bp) = 5ap. In this 
case the matrix V (often called Gram matrix) is equal to the 
identity matrix, i.e., V = I, and it follows that 

x + s = WV s = W s. (7) 

Note that our simplifying assumption of orthogonal basis 
functions can be easily relaxed, leading to slightly modified 
W —WV and C —?■ CV matrices (below). 

The numerical evaluation of Galerkin inner products in 
Eq. 6 involves multidimensional integrals over the support 
of both basis and testing functions. This integration can be 
quite cumbersome due to singularities (when the support of 
the basis and the testing functions overlap) and the highly di¬ 
mensional aspect of the problem. However, previous work*®^ 
demonstrated that these challenging volumetric integrals can 
be reduced to surface integrals (of lower singularity), allow¬ 
ing us to benefit from decades of work dedicated to the ac¬ 
curate and efficient evaluation of the associated surface inte¬ 
grals. Here, we make use of the free-software DEMCEM'^^ 
and DIRECTFN,'^^ which leverage the techniques described 
in Refs. 162 and 163 . Furthermore, MoM JM-VIE formu¬ 
lations with local basis/testing functions typically result in 
very large linear systems, which can be solved with iterative 
algorithms for non-symmetric dense systems. In each itera¬ 
tion, the associated matrix-vector products take 0{N‘^) time. 
Moreover, it is practically impossible to explicitly store the 
(dense) matrix W~^ requiring 0{N'^) memory. In fact, there 
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where the superscripts denote blocks associated with the vari¬ 
ous objects, with diagonal components corresponding to self¬ 
interactions and off-diagonal blocks involving interactions be¬ 
tween different objects. Finally, we define the projection. 


pp 

a,p 


1, if a = P = p 
0, otherwise, 


( 10 ) 


which selects specific blocks of vectors x^ = P^x or diagonal 
blocks of matrices Ap = PPAPP corresponding to object p. 


B. Power transfer 

We now derive a compact matrix-trace formula for the 
computation of the ensemble-averaged flux into body Bp (or 
equivalently the absorbed power) due to fluctuating current 
sources in body Bq, integrated over all possible positions and 
orientations. The first step consists of the evaluation of the 
flux from Bp due to a single dipole source a immersed in Bq, 
which we denote as . Direct application of Poynting’s 
theorem implies that the flux on the objects is given by: 

<f>yp = ^Re[ d^^e-4> ( 11 ) 

JBp 

which amounts to the work done by the total held on the po¬ 
larization currents in Bq. Expressing the induced currents and 
fields in the basis of JM-VIE currents and using the relation 
(^ = r * (^ -f cr) yields the following discrete approximation 
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(see Ref. 155 for a complete analysis): 

= i RexP*(l)P = i Rex*PP(l) 

= ^Re{x + s‘^ypPG{x + s‘^) 

= i(x + s'^)* sym (PPG')(x + x'*) (12) 

= i {WP’^s)* sym {PPG) {WP‘^s) 

= ixr [{ss*)(WPy* sym (pPG){WPy] 


where sym G = denotes the Hermitian part of G. It 

is then straightforward to obtain the ensemble-averaged flux 
which yields: 

= ^Tr [{ss*){WPy* sym{PPG){WPy] 

? (13) 

= -Tr [P'^GP'^W* sym {PPG)W] 

where C = (ss*) is a current-current correlation matrix 
that captures a statistical, ensemble average over sources, 
described in more detail in Sec. HE. Defining the matrix 
— pq(jpq^ which is simply a projection of the correla¬ 
tion matrix unto the space of basis functions in q, we find that 
the ensemble-averaged flux is given by: 


2 


C^W* sym {PPG)W 


(14) 


about some origin Xq can be obtained by integrating the dif¬ 
ferential torque dr = (x — Xq) x dF on a volume element. 

Expressing the induced currents and fields in the basis of 
JM-VIE currents and following a similar procedure as that of 
Sec. IIB, one finds that the ensemble-averaged force on the 
object can be written in the compact and convenient form: 


jrg->-p _ _-pj. 

2a; 


C'^W* asym (P^G^) W , 


(16) 


where in this case and in contrast to power transfer, the rele¬ 
vant quantity is the matrix representation G^ of the gradient 
of the Green’s function operator G, whose matrix elements 
^a ,/3 = (ba,Vr * bjs). Also, asymG = denotes 

the skew-Hermitian part of G. The torque on the object can 
be obtained similarly by computing angular derivatives of G. 
It turns out that the calculation of these matrix elements re¬ 
quires evaluating multidimensional integrals whose singular¬ 
ities are more severe than those of G. A key distinction be¬ 
tween fluctuation-induced transfers of power and momentum 
is that, in the latter case, one finds nonzero fluctuation-induced 
forces and torques between bodies even at thermal equilibrium 
and even at zero temperature; these are just the usual equilib¬ 
rium Casimir forces.Equation 16, which computes only the 
non-equilibrium contribution to the force, must generally be 
augmented by these equilibrium contributions to yield the to¬ 
tal force. Connections between Eq. 16 and expressions for 
equilibrium forces, along with techniques for evaluating the 
above-mentioned integrals and results of VIE computations 
of non-equilibrium Casimir forces and torques are addressed 
in subsequent work.*®^ 


C. Momentum transfer 

In addition to carrying energy, the radiation emitted by fluc¬ 
tuating sources also carries linear and angular momentum, 
which can also be described using similar expressions. The 
starting point consists of the evaluation of the force (or torque) 
imparted on an object Bp due to a single dipole source im¬ 
mersed in Bq. Although electromagnetic forces are often com¬ 
puted via surface-integrals of the Maxwell stress tensor, it is 
also possible and in our case more convenient to express the 
force as a volume integral by considering the Lorentz force 
acting on the internal currents ^ induced on In particu¬ 

lar, the force on the object is given by: 

= (15) 

JBp 

where V denotes the usual partial derivative with respect to 
infinitesimal displacements. The derivation of the above ex¬ 
pression follows from application of the time-average Lorentz 
force dF = ^Re (p*E -f J* x B)c?^x on the electric charge 
and current densities (p, J) in an infinitesimal volume element 
d^x, together with a similar expression for the force on the 
magnetic sources. Integrating over the volume of the body and 
employing Stokes’ theorem along with Maxwell’s equations 
immediately yields Eq. 15. In a similar fashion, the torque 


D. Far-fleld radiation intensity 

In addition to power and momentum transfer, another use¬ 
ful quantity is the far-fleld radiation intensity of our system, 
which can also be expressed as a simple trace formula. The re¬ 
sult which follows trivially from Eq. 13, is that the ensemble- 
averaged flux radiated by an isolated body Bq to the back¬ 
ground medium is given by: 

=-^Tr [GW* sym GW] (17) 

where the minus sign corresponds to the direction of the power 
flux and stems from Poynting’s theorem. However, in addition 
to the overall radiation, it is also useful to obtain the radiation 
intensity over specific directions, or equivalently the power ra¬ 
diated per solid angle. The angle-resolved radiation intensity 
fj-om a single source a immersed in Bq can be obtained 
by expressing the radiation held at infinity Eoo (where only 
far held contributions remain) in terms of the free and bound 
current sources, as follows: 

1,2 y i-2 7 

(18) 

where k is the wavenumber and Z = y/po/eo is the wave 
impedance, both in vacuum. Also, r^(x,y) is the 3x6 
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Green’s tensor of the ambient medium which maps currents 
to far-held electric helds, and Q is a 3 x 3 transformation ten¬ 
sor that maps vectors from Cartesian to spherical coordinates 
and projects their radial component to zero.Given the so¬ 
lution of the VIE scattering problem and following the same 
procedure described above, it is straightforward to write the 
radiation intensity as a matrix-trace formula of the form: 

Li (19) 

= 2^Tr [(ss*)(IE)*(G^*G^)(IV)] 

where the matrix G^ is the discretized form of the operator 
obtained in a similar fashion as G. Ensemble averaging 
over all sources, we hnd that the hnal formula for the angle- 
resolved radiation intensity is given by: 


( 20 ) 


Equation 20 can be integrated over all solid angles fl to yield 
the total radiation rate = J dfl which as ex¬ 

pected agrees with results obtained by direct application of 
Eq. 17, as discussed in Sec. III. 


E. Current-current correlation matrices 

The formulas above are very general in that they apply to 
many different kinds of huctuation processes, the physical 
properties and origins of which are embedded in the correla¬ 
tion matrices G = (ss*), involving ensemble averages over all 
sources a and polarizations throughout the bodies. In partic¬ 
ular, the matrix elements of the correlation matrices describe 
interactions among basis functions and are given by: 

= (saSp) = J J (fiyid^yb*^{x){a{x)a*{y))bfs{y) 

( 21 ) 

which follows trivially from the orthogonality property of our 
basis functions and the fact that cr(x) = ^abaip^)- Al¬ 
though in general the calculation of each matrix element in¬ 
volves volume-volume integrals against pairs of basis func¬ 
tions, current fluctuations are temporally and spatially uncor¬ 
related in local media^’^’^° and are described by: 

{a,{yi,uj)a*{y,uj)) = Jy (x, w)(5(x - y) (22) 

where the subscripts denote polarization degrees of freedom 
and ^7 > 0 is a position-dependent spectral tensor whose form 
depends on the physical origins of the fluctuations. It follows 
that G is Hermitian and positive-semidefinite and thus admits 
a Cholesky factorization G = LqLq, which we exploit in 
Sec. Ill to demonstrate that our radiation, power, and momen¬ 
tum formulas are susceptible to fast-trace calculations. 


When the sources of fluctuations involve only quantum and 
thermal vibrations (heat), the correlation function J is de¬ 
termined by thermodynamic considerations such as the well- 
known relating current fluctuations to dissipation 

in materials. Without loss of generality, the spectral function 
is given by:'™ 


4 

J)j (x,w) = - Imxy (x, w)0(x,w), (23) 

TT 

where the Imy tensor describes losses in the medium and 
0(x, w) = huj— 1) is the Planck distribution, 
or the average energy of an oscillator having local tempera¬ 
ture T(x). Equation 23 in conjunction with the power transfer 
and radiation formulas above are exploited below to evaluate 
thermal radiation and heat transfer between inhomogeneous 
bodies with spatially varying temperature and dielectric prop¬ 
erties, and also in an upcoming paper that focuses on non¬ 
equilibrium Casimir forces.'®^ 

In situations involving active media driven by external 
pumps, the properties of the fluctuating currents and hence J 
depend on the details of the input drive along with the physi¬ 
cal emission mechanisms. Eor a broad range of processes, the 
spectral function can be written in the simple form: 


i7ij(x, w) — yiiic(x)yemm,ij (X) w)) (24) 

where Xinc describes the response of the medium due to the 
pump and Xemm describes the emission spectrum of the ex¬ 
cited medium, which depends on the distribution of active 
molecules in the medium and on complicated electronic tran¬ 
sitions mediated by the pump as well as quantum/thermal pro¬ 
cesses.^ In the particular example of one-photon fluorescence 
from a medium (with high quantum yield) excited by inci¬ 
dent light, the pump spectrum is proportional to the locally 
absorbed power and hence can be computed by direct appli¬ 
cation of the VIE power formulas. Such a relationship in con¬ 
junction with Eq. 20 is exploited below to compute the fluo¬ 
rescence spectrum of an irradiated sphere. A similar depen¬ 
dence on the local held intensity arises in the case of Raman 
scattering, except that Xinc is proportional to the Raman polar¬ 
izability tensor rather than the susceptibility of the medium.^ 
In the case of spontaneous emission from a gain medium, the 
emission spectrum is determined by spatially dependent effec¬ 
tive permittivity and temperature profiles determined by the 
driven steady-state atomic populations of the medium, both of 
which can be obtained by application of steady-state ab-initio 
laser theory (SALT).^’'7i Similar descriptions apply in more 
complicated systems, including fluorophores with low quan¬ 
tum yields or active media subject to highly nonlinear (e.g. 
two-photon) processes. 


III. FAST TRACE COMPUTATIONS 


The matrix-trace formulas derived in the previous sections 
require products of inverses of the JM-VIE matrix W with 
dense matrices sym {P^G), asym (P^VG), and G^*G^. As 
mentioned above, due to their large size and correspondingly 
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severe CPU and memory limitations, it is practically impos¬ 
sible to form explicitly either the Green matrix or its inverse. 
There are however fast FFT-based procedures for evaluating 
matrix-vector products of the JM-VIE system matrix and the 
Green matrix.Here we describe a framework based on iter¬ 
ative methods for the fast computation of the associated trace 
formulas above. 

We begin with the matrix-trace formula <1)5“*^^ in the pres¬ 
ence of n bodies (including Bp and Bq), which after some 
algebraic manipulations can be written as follows (ignoring 
pre-factors): 


= Tr [C««fFP«*(symGPP)lUP«] 


-b ^ Tr [G«« sym (25) 

m—l 

m^p 


where is the qq block of the matrix C. Due to the dif¬ 
ferent characteristics of and we need to address 

them separately. As discussed in Sec. HE, the matrix C‘^‘^ can 
be assumed to be Hermitian and positive semidefinite, hence 
it admits a Cholesky factorization, = Lqii L’^qq ■ In addi¬ 
tion, symG^’^ is a Hermitian, negative semidefinite matrix 
and it also admits a low-rank approximation since it is asso¬ 
ciated with the smooth, imaginary part of the Green’s func¬ 
tions. Hence, it can be approximated to any desired accuracy 
by a truncated singular value decomposition (SVD) factoriza¬ 
tion, sym GPP ss -UppSppUpp*, where S'pp G with 

r "C iV. The norm of the error in the aforementioned trunca¬ 
tion is bounded by the norm of the vector of discarded singu¬ 
lar values. The classical SVD algorithm requires the complete 
matrix, hence we resort here to a class of modern randomized 
matrix approximation techniques, and more specifically to the 
randomized SVD method (rSVD).'^^’*^"^ rSVD is effective for 
matrices with fast drop of the singular values and it requires 
only a fast matrix-vector procedure, which we have developed 
as described above. The matrix with the singular values can be 
further decomposed so that Spp = LsppLgpp. Einally, it fol¬ 
lows that the self-term in Eq. 25 can be written as the square 
of a Erobenius norm, 

S<i^p = -TT[LcppL*cqq{WP‘^*UPP)LsPpL*spp{UPP*WP‘>)] 
= -\\L*Cqq{WP‘>*UPP)LsPp\\l 

(26) 

Eor the most time consuming part of the norm, we need to 
solve the adjoint JM-VIE system r times (for each of the lead¬ 
ing singular vectors of symGPP). Note however that we can 
solve for each vector of IJPp independently and thus the entire 
procedure is embarrassingly parallelizable. Also, L^qq and 
Ls PP are either sparse or diagonal, while WP‘^*Upp is a “tall- 
and-skinny” matrix (the number of columns is much smaller 
than the number of rows) and hence the matrix product ap¬ 
pearing in the norm can be computed efficiently. 

The trace formula for C^^p is not symmetrical and there¬ 
fore cannot be reduced to a norm. In this case, one can ex¬ 
ploit the fact that Gp™ admits a low-rank approximation due 


to the smoothing properties of the Green’s function for dis¬ 
joint objects. The final dimensions of the low-rank approx¬ 
imation of C‘^^P (for a prescribed accuracy) depend on the 
electric distance between objects p and i.e., GP'^ ~ 

upmgpmypm*^ where 5'P™ G with I < N. The fi¬ 

nal formula for G^^p after the Cholesky factorization of the 
singular values matrix (Sp^) is given by 

n 

G'J^P = Re ^ Tr [Xupr^X^^pJ^ ^21) 

m—l 

m^p 


where 

Xupm = L*Cqq{WP'^*UP"^)LsprP 
Xvpm = L^„(lU'"'?*UP™)Tsp™. 

Both Xjjp^ and Xypm are “tall-and-skinny”, and we can not 
compute the trace by forming explicitly their product, due to 
memory limitations. Alternatively, we can use the standard 
vectorization of a matrix vec(), which converts the matrix 
into a column vector, together with the identity, Tr [XF*] = 
vec(X)^ • vec(F), and write Eq. 27 in the following compu¬ 
tationally friendly form: 

n 

(jq^p — p(g ^ vec(X(7pm)^ • vec(Xvpm). (2g^ 

m—l 

m^p 


The overall computational complexity for the evaluation of 
G<i-^p consists of a single run of the Randomized-SVD for a 
non-symmetric matrix,and 2 x I solves of the adjoint JM- 
VIE system. In the case of the matrix-trace formulas for the 
force and the torque, the procedure is similar with the one 
described above. The only difference stems from the replace¬ 
ment of G with G^ and sym with asym. 

Einally, the case of far-held radiation is somewhat simpler. 
According to Eq. 20, we just need to solve 2 times the adjoint 
JM-VIE system, since G^ G Hence, the radiation 

intensity for a specihc direction or solid angle H, is given by 
the following square of the Erobenius norm: 

(29) 

This is a very useful formula, especially when directional 
information of the radiated power is of interest. In addi¬ 
tion, the total radiated power can be evaluated by integrat¬ 
ing Eq. 29 over all solid angles, as mentioned in Sec. IID, 
which would amount to employing a numerical integration 
scheme over the unit sphere (e.g. Lebedev quadrature'^®). 
Alternatively, one could exploit Eq. 17 and the associated 
norm \\Lq{W*U)Ls\\p to compute the total radiated power 
from an isolated body. The latter is expected to be more ef- 
hcient for total-radiation computations with prescribed accu¬ 
racy, controlled by the SVD factorization of the Green matrix, 
in which case the minimum number of JM-VIE solves needed 
for a prescribed accuracy is estimated in advance. In contrast, 
the former approach is based on adaptive quadrature schemes 
where the accuracy is controlled by the comparison of results 
between different orders of integration, with no a priori con¬ 
trol. 
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Figure 2. Flux spectrum <1 ?(cj) of a cube of edge-length 2R held at 
temperature T, normalized by the corresponding black-body spec¬ 
trum T), for different (a) discretization 

mesh densities and (b) rSVD truncation tolerances. 



IV. VALIDATION AND APPLICATIONS 

In this section, we apply the FVC method to obtain new re¬ 
sults in complex geometries. To begin with, we show that 
the Green matrices appearing in our trace formulas admit 
low-rank decompositions (as discussed in Sec. Ill) by com¬ 
puting their ranks to within some tolerance in a representa¬ 
tive structure involving two vacuum-separated, homogeneous 
cubes. We validate the FVC method by checking its pre¬ 
dictions against known results of thermal radiation and near- 
held heat transfer between homogeneous bodies, including 
spheres, cubes, and ellipsoids, obtained using a boundary- 
element implementation of our recent FSC formulation.^® 
We show that when subject to temperature gradients or con¬ 
tinuously varying permittivities, complex bodies can exhibit 
highly modihed thermal radiation and heat transfer spectra, 
leading to directional emission at selective wavelengths. Fi¬ 
nally, we demonstrate that the same formalism can be ex¬ 
ploited to study luminescence from excited media by com¬ 
puting the Huorescence spectrum of a sphere irradiated by 
monochromatic incident light. We show that the impact of the 
resulting inhomogeneous current Huctuations cannot be eas¬ 
ily obtained by exploiting simple homogenization or effective- 
medium approximations. For convenience and simplicity, we 
consider dielectric media with no material dispersion (con¬ 
stant Ree « 12 and large dissipation Ime Ri 1), though our 
approach is general in that it can readily handle other kinds of 
materials such as metals with Re e < 0 and even gain media. 


A. Low-rank approximations 

Low-rank approximations of the associated (free-space) 
Green matrices are instrumental to the practical and effi¬ 
cient evaluation of our trace formulas. In this section, we 
present some representative results obtained from computing 
the ranks of both sym and to within some toler¬ 

ance, for the particular problem of two vacuum-separated, ho¬ 
mogeneous cubes of edge-length L = 2R and separated by a 
surface-surface distance d, shown schematically in Fig. 5. 


Table I. Ranks of symG^^ for various frequencies and toler¬ 

ances (tol) in truncated SVD. The ranks correspond to the case of a 
cube of edge-length 2R. In addition, results for a sphere of radius R 
are included in brackets. 


tol 

C 

le-l 

le-2 

le-3 

le-4 le-® 

le-® 

0.01 

4(4) 

4(4) 

4(4) 

4 (4) 7 (7) 

12 (12) 

0.1 

4(4) 

4(4) 

7(7) 

12(12) 12(12) 

14(12) 

1.0 

12(7) 

14 (12) 24 (24) 40 (24) 40 (40) 

60 (40) 

2.0 

18 (12) 

37(24) 51 (40) 

65 (60) 84(60) 

109 (84) 


Table II. Ranks of for various distances (d) and tolerances (tol) 
in truncated SVD. The ranks correspond to the case of two cubes 
of edge length L = 2R and frequency — 1. Each cube is 
discretized into N = 40^ voxels, resulting in 3N total degrees of 
freedom, i.e., #DOFS = 3N. 


tol 

le-l 

le-2 

le-3 

le-^ 

le-S 

le-® 

0.001 

4075 

4853 

5253 

6352 

7240 

8481 

0.01 

992 

2611 

3934 

4800 

5832 

6894 

0.1 

50 

196 

447 

804 

1268 

1849 

1.0 

6 

14 

27 

42 

66 

89 

10.0 

4 

7 

9 

14 

19 

23 


Table I shows the singular values of sym correspond¬ 
ing to one of the two cubes, as a function of the normalized 
frequency coRjc and tolerance tol; that is, we obtain the sin¬ 
gular values that produce SVD factorizations bounded in norm 
by the tolerance tol, also known as a truncated SVD. Since the 
associated matrix is very large and our trace formulations can 
be cast in terms of fast matrix-vector products, our calcula¬ 
tions exploit the rS VD method recently developed for big-data 
problems. (Note that results for the second cube, involving 
sym would be identical since both cubes have equal sizes 
and number of unknowns.) Our results reveal at least two im¬ 
portant features; First, the ranks scale linearly with ui at large 
frequencies, and sub-linearly (roughly constant) at small fre¬ 
quencies. Additional numerical experiments (not shown) con¬ 
firm that the effect of mesh density on the ranks is negligible, 
yet another manifestation of the favourable convergence prop¬ 
erties of the JM-VIE formulation.^® This also suggests a strat¬ 
egy for obtaining the finite rank of sym G^p with prescribed 
accuracy; we begin by computing the rank of the operator for 
a prescribed accuracy by using a coarse mesh and then run a 
fixed-rank rSVD algorithm with finer mesh. Finally, Fig. 2 
illustrates the rate of convergence of the radiation spectrum 
<h(a;) from an isolated cube at a fixed temperature T with 
respect to different (a) discretization mesh densities and (b) 
truncation tolerance, normalized to the spectrum of a corre¬ 
sponding black body <1'bb(w) = ^^(a;/c)^0(u;,T), where 
A denotes the surface area of the cube. 

The situation changes in the case of the “coupling” Green 
matrix G®^, which encodes interactions between objects. Ta- 
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ble II shows the significant singular values associated with the 
coupling matrix of the same cube-cube geometry at a fixed 
frequency uj and for various separations d, obtained by lever¬ 
aging the rSVD technique. As expected, the singular values 
increase as d decreases, a consequence of the power-law drop¬ 
off of the Green’s function with separation in the near field. It 
follows that the computation complexity of the trace formu¬ 
las increases as the two bodies come close together. (Note 
that, as described in Sec. Ill, our trace formulas for power and 
momentum transfer require us to solve two VIE systems for 
every corresponding eigenvector, but fortunately each system 
can be solved independently and the overall process is em¬ 
barrassingly parallelizable.) Nevertheless, we find that 
remains very low rank even for relatively close separations 
d/L Ki 0.1, below which constraints on the resolution make 
the FVC approach less practical. However, it is precisely at 
such small separations that approximate methods such as the 
proximity approximation become accurate.^* 


B. Thermal radiation and heat transfer 

We begin by validating our FVC approach by checking 
its predictions of thermal radiation from homogeneous bod¬ 
ies against results obtained using our recently developed FSC 
formulation,^®’^^ which is well-suited for handling piece-wise 
constant structures and fluctuations statistics. Figure 3(a) 
shows the flux spectra $(aj) of multiple objects (of uniform 
temperature T and permittivity e = 12 +i, including a sphere 
of radius R (blue line), a cube of edge-length 2R (green line), 
and an prolate ellipsoid of long semi-axis R and short semi¬ 
axes (red line). Note that in each case is normalized 
to the corresponding flux from a black body. As shown, there 
is excellent agreement between the FVC (solid lines) and FSC 
(circles) predictions, both of which illustrate the expected ra¬ 
diation enhancement at geometric resonances. 

The FVC method can also handle more complex structures, 
including inhomogeneous bodies with spatially varying per¬ 
mittivities. In particular. Fig. 3(b) shows <I>(w) for the same 
geometries of Fig. 3(a) but for objects with linearly varying 
permittivity profiles e{z) = t-R -b (e_R — with 

= 2-|-fandei^ = 12-|-f (solid lines) and axes chosen to lie 
at the geometric center of each object. Compared to the spec¬ 
trum of the homogeneous bodies of Fig. 3(a), one finds that 
the resonances are shifted to larger frequencies and their peak 
amplitudes are significantly smaller, a consequence of the de¬ 
creased effective permittivity of each object. For comparison, 
we also show <I>eff(w) (dashed lines) from corresponding ho¬ 
mogeneous objects with effective permittivities, 

Ceff = [ d^xejx), (30) 

k Jv 

corresponding to uniform Ceff = 7 + i. Our calculations re¬ 
veal that in the illustrated frequency range and for our choice 
of dielectric profiles, the homogeneous approximation is qual¬ 
itatively accurate to within 10%. On the other hand, employ¬ 
ing Eq. 20 to compute the angular radiation patterns at se¬ 
lected frequencies, shown as insets in Fig. 3, reveals signif- 



Figure 3. Flux spectrum <l?(a;) normalized by the corresponding 
black-body spectrum <I>BB(t<j) = j^(cu/c)^0(cj, T) of different 
bodies of surface area A held at temperature T = 1000 K, including 
a sphere of radius R (blue lines), cube of edge-length 2R (green line), 
and ellipsoid of long semi-axis R and short semi-axis A (j-gd line). 
The objects have either (a) uniform permittivities e = 12 -|- i or (b) 
spatially varying e(z) = e_i?-b (er? — e_fl)with erj = 12-bi 
and = 2 -b i. For comparison, we also plot the radiation spec¬ 
trum <l>eff (dashed lines) of corresponding bodies with homogeneous 
effective permittivities Ceff = 7 + i. The insets depict the angular 
distribution of far-field radiation U (H), normalized by the maximum 
intensity over all directions maxn U, at selected frequencies. 


icant changes, e.g. significantly larger directional emission, 
that cannot be captured by the effective-medium approxima¬ 
tion. In particular, the radiation patterns of the inhomoge¬ 
neous objects break z mirror symmetry. For example, the flux 
from the cube at w Ri 0.65i?/c is slightly larger in the —z 
than in the +z direction, a situation that is reversed at larger 
u! Ri 0.9i?/c (see insets). Generally, the transition frequency 
of the favored radiation direction depends on the geometry; 
for instance, even at a frequency as high as w « 1.5i?/c, the 
ellipsoid continues to radiate more along the —z direction. 

More pronounced changes arise when objects are subject to 
spatial temperature gradients. Figure 4 shows from ho¬ 
mogeneous (e = 12-bf) ellipsoids subject to either (a) radially 
varying T(r) = Tq + (Tr — Tq)-^ or (b) z-varying tempera- 
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Figure 4. Flux spectrum 4?(cu) of various bodies normalized by 
the corresponding predictions of a simple approximation $eff, de¬ 
fined in Eq. 31, including (a) sphere of radius R and radially vary¬ 
ing temperature profile T{r) = To + {Tr — Tb)-^ for both To = 
0, Tr = 1000 K (blue line) and To = 1000, Tr = 0 (green 
line), and (b) sphere of radius R (blue line) or ellipsoids with short 
semi-axes y and long semi-axis R along the z (green line) or x 
(red line) directions, subject to vertically varying temperature pro¬ 
files T{z) — T-l + (Tl — where L denotes the 

2 -dimension of the corresponding body. In all cases, objects have 
uniform permittivity e = 12 -f i and are subject to temperature gra¬ 
dients T-l = 0 and Tl = 1000 K. The insets in (a) show the local 
density of states along a cross-section of the sphere at different fre¬ 
quencies while those in (b) show the angular distribution of far-field 
radiation U(yi) normalized by maxn U. 


ture profiles (see caption). In both cases, <!> is normalized by 
the flux $eff obtained from a naive approximation in which 
the temperature variations are removed in favor of a uniform 
effective temperature determined by a simple average of 
the Planck distribution over the volume V of the bodies, 

Q{uj,Tee) = ^ (31) 

Such a simple approximation obviates the need for exact cal¬ 
culations that explicitly incorporate inhomogeneities, but is 
clearly inadequate for wavelength-scale objects. Specifically, 
Fig. 4(a) shows from spheres with radially varying tem- 



Figure 5. Heat-transfer spectrum $(aj), normalized by the corre¬ 
sponding black-body spectrum 4 ?bb(uj) = ■^^(oj/c)^Q(cj,T), be¬ 
tween two cubes of edge-length 2R and temperature T = 1000 K 
separated by surface-surface distance d = R. The cubes are as¬ 
sumed to have either uniform permittivities e = 2 + i (red dashed 
line), e = teff = 7 + i (black dashed line), or e = 12 4- i (blue 
dashed line), or vertically varying permittivities e( 2 i) = e-R + {€R — 
s-r) defined with respect to the local axis xi ,2 at the center of 

each cube (shown on the inset), chosen so that the system has mirror 
symmetry about the x-y plane intersecting the origin O. The gradi¬ 
ents are either increasing (black solid line) or decreasing (green solid 
line) toward or away from the center, corresponding to the choice of 
CH.-i? = {12 + i,2 + i} or eR ^ e-R, respectively. 


peratures, illustrating that beyond the sub-wavelength regime 
oj ^ R/c and depending on the choice of Tq and Tr, $ can 
be many times larger or smaller than that predicted by Eq. 31. 
The failure of this naive approximation is especially appar¬ 
ent near resonances, where the coupling of fluctuating sources 
(dipoles) to far-held radiation (the local density of states) is 
highly position-dependent. The insets of Fig. 4(a) show cross- 
sections of the spatially varying fiux contribution from dipoles 
in the interior of the sphere at two relatively close frequen¬ 
cies. At ujR/c ~ 1.1, we find that dipoles closer to the center 
can couple more efficiently to far-field radiation than those 
near the edges, causing Eq. 31 to underestimate the flux by 
<l>/<l>eff 3 in the case Tq = 0, Tr = 1000 K (green line) 
and to overestimate it by <l>/<l>eff w 0.8 when Tq = 1000 K, 
Tr = 0 (blue line). The converse is true at ujR/c si 0.85, in 
which case their coupling to radiation is largest at the center 
and edges of the sphere. Similar effects arise in situations in¬ 
volving z-varying temperature profiles, explored in Fig. 4(b) 
for either spheres (blue line) or ellipsoids with either their 
long-axes (green line) or short-axes (red line) aligned with the 
z direction. For instance, ellipsoids can exhibit highly direc¬ 
tional emission (almost a factor of 3 times larger) along the 
direction of increasing temperature. 

In addition to far-field radiation, the FVC method can be 
employed to obtain radiative transfer between objects. Fig¬ 
ure 5 shows the heat-transfer spectrum $(a;) (computed via 
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Eq. 16) normalized by <i>BB(w) (same as above), between 
two vacuum-separated cubes of edge-length 2i? and surface- 
surface separation d = R, of either uniform (dashed lines) or 
vertically varying (solid lines) permittivities. We consider di¬ 
electric prohles of the form e{zi) = c-r + {cr — 
dehned with respect to the local axis located at the center 
of each cube xi 2 , chosen so that the entire system has mir¬ 
ror symmetry about the origin (see inset). We consider two 
different prohles, C-r^r = {2 + i,12 + i] (black line) or 
tR -H- e-R (green line), corresponding to increasing gradi¬ 
ents toward or away from the origin. For comparison, we 
also plot the transfer between cubes of uniform permittivities 
e = 2-|-i (red dashed line), e = 12-|-i (green dashed line), and 
e = Cefi = Y corresponding to the minimum, 

maximum, or average of the spatially varying permittivities, 
respectively. As shown, depending on the wavelength regime 
(near versus far held) inhomogeneities can have a different ef¬ 
fect on the heat transer. For instance, at low wi?/c 1 where 
near-held effects prevail, homogeneous bodies with smaller 
dielectric constants tend to transfer more heat—the same de¬ 
pendence is observed for planar objects separated by vacuum, 
where the near-held contribution ~ Not surpris¬ 

ingly, because nearby regions tend to contribute more than 
far-away regions, one observes that despite having the same 
average permittivities CeS (dashed blue line), the transfer is 
sensitive to the local dielectric variation, exhibiting larger en¬ 
hancement in the case where the permittivity is increasing to¬ 
ward (green solid line) rather than away (black solid line) from 
the origin. At larger ujRjc > 0.5 where far-held effects be¬ 
gin to dominate, one observes the opposite behavior, in which 
case the largest transfer is obtained for decreasing permittivi¬ 
ties toward the origin. Essentially, as illustrated in Fig. 3(b), at 
sufficiently large wavelengths, bodies with dielectric gradients 
tend to radiate along the direction of increasing permittivity. 


C. Fluorescence 

We now consider application of the FVC formulas to the 
calculation of huorescence. A typical huorescence setup con¬ 
sists of an incident wave impinging on a huorescent body, 
leading to the absorption and subsequent re-emission of light 
by molecules inside the body.^ Both of these effects are cap¬ 
tured by the current-current correlation matrix described in 
Sec. HE, which encodes the spectral properties of the huctu- 
ations. In the particular problem of one-photon huorescence 
induced by an incident monochromatic wave at a given fre¬ 
quency Wine, the spectral function J (x, w) has the form given 
in Eq. 24, with the excitation spectrum given by the locally 
absorbed power, 

Xinc (x) cx WincImxlE (x,Wi„e)P, (32) 

and Xemm(x, w) denoting the Huorescence spectrum of the 
bulk medium, usually a relatively broad Lorentzian lineshape 
centered near the material’s absorption resonance. (Note that 
Xinc = 0 in the absence of a huorescent medium.) A well- 
known approach to enhance huorescence involves design- 


X-Y X-Z Y-Z 
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Figure 6. Far-held fluorescence spectrum <l?(a;) (in arbitrary units) 
of a homogeneous and non-dispersive dielectric sphere of radius R 
and permittivity t = 12 -\- i excited by an x-polarized planewave 
propagating along the z direction with frequency uJi^cR/c = 1.58. 
The absorbed power Xinc (x) inside the sphere, obtained by solving a 
single scattering problem as described in Ref. 25, is shown in the top 
contour plots along three sphere cross-sections. $ is computed ex¬ 
actly (blue line) or via a homogeneous approximation <l?eff in which 
the absorbed power is taken to be uniformly distributed inside the 
sphere and given by Xeff = d^xxinc(x) (red line). The ratio of 
the two is plotted as the black dashed line on the right axis. The insets 
depict the angular distribution of fluorescence emission, normalized 
by the maximum intensity over all directions, at selected frequencies. 


ing bodies to have strong resonances at Winc, leading to in¬ 
creased absorption.^ For bodies designed to have additional 
resonances within the fluorescence bandwidth, determined by 
Xemm, there is an additional source of enhancement arising 
from the increased local density of states, or increased cou¬ 
pling of dipole emitters to far-held radiation. Inhomogeneities 
arise due to the fact that xinc and the local density of states are 
both highly spatially non-uniform near resonances. 

Figure 6 shows the fluorescence emission <l)(a;) from a 
sphere of radius R and uniform permittivity e = 12 -F i, ir¬ 
radiated by an x-polarized, z-traveling incident wave of fre¬ 
quency uJincR/c ~ 1.58, chosen to coincide with one of its 
resonances. For simplicity, we assume a non-dispersive and 
uniformly distributed fluorescent medium with Xemm = 1, al¬ 
though as noted above our formalism can just as easily handle 
spatially varying distributions. The flrst step in computing the 
fluorescence emission is to obtain the locally absorbed power 
within the sphere Xinc(x), which boils down to the calcula¬ 
tion of a single and far simpler scattering problem exploiting 
Eq. 12, as described in Ref. 155. Along with $ (blue line). 
Fig. 6 shows Xinc along three different cross-sections inter¬ 
secting the center of the sphere (top contour plots), illustrat¬ 
ing the highly non-uniform spatial pattern of current fluctu- 
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ations. Also shown is the spectrum $eff obtained by appli¬ 
cation of a homogeneous approximation (red line) where the 
absorbed power is averaged over the volume of the sphere to 
yield a uniform, effective Xeff = fy Xinc(x), along with 
the corresponding ratio <l’/$efr (black line). As before, such 
approximations yield accurate results in the sub-wavelength 
regime but break down at larger frequencies. For instance, at 
ujR/c 1 we find that $/3>eff ~ 1-5. More importantly, 
the approximation fails to capture the angular distribution of 
radiation (insets): both the direction of largest fluorescence 
and overall emission pattern change drastically as the emis¬ 
sion frequency increases from uiR/c ~ 1.1 to loRjc « 1.3. 

V. CONCLUDING REMARKS 

Our FVC formulation of electromagnetic fluctuations en¬ 
ables accurate calculations of wide-ranging incandescence 
(e.g. thermal radiation, dispersion forces, heat transfer) and 
luminescence (e.g. spontaneous emission, fluorescence, Ra¬ 
man scattering) phenomena in arbitrary geometries. Simi¬ 
lar to recently proposed scattering-matrix and surface-integral 
equation formulations of radiative heat transfer, the resulting 
quantities are obtained via traces of matrices involving inter¬ 
actions among basis functions; however, because the JM-VIE 
“scattering” unknowns are volume currents rather than prop¬ 
agating waves or surface currents, the formalism is applica¬ 


ble to a broader set of problems. For example, as demon¬ 
strated here, our approach captures phenomena associated 
with the presence material inhomogeneities, such as spatially 
varying temperature gradients and dielectric properties within 
bodies. In future work, we plan to exploit the FVC ap¬ 
proach to demonstrate predictions of highly directional radia¬ 
tion from inhomogeneous structures subject to thermal gra¬ 
dients,non-equilibrium Casimir torques on chiral parti¬ 
cles,'®^ and enhanced directional emission from parity-time 
symmetric (gain) media.Furthermore, although our calcu¬ 
lations focused on geometries involving compact bodies, the 
same power and momentum formulas derived above apply to 
geometries involving extended bodies, the subject of future 
work. 
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